---
title: "ISQ-R&R"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

pacman::p_load(foreign,
               tidyverse,
               texreg, 
               stargazer,
               ggrepel,
               kableExtra, 
               corrplot, 
               ggcorrplot,
               Hmisc, 
               styler, #highlight a code chunk and press CNTRL+SHIFT+A for formating 
               plm,
               biglm,
               lfe,
               haven, 
               knitr)

```


#Loading the data
```{r}
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#loading data
naturaldis <- read_csv("natdis_june2018.csv")

```


### Table 2 ###
```{r testt, }

naturaldis2 <- read_csv("natdis_august2020.csv")

depth <- felm(data=naturaldis2, ln_exports ~  deeppta + naturaldis_deeppta | ID1 + ID2 + ID_dyad)
flex <- felm(data=naturaldis2, ln_exports ~  flexpta + naturaldis_flexpta | ID1 + ID2 + ID_dyad)
enf <- felm(data=naturaldis2, ln_exports ~  enfpta + naturaldis_enfpta | ID1 + ID2 + ID_dyad)
contdepth <- felm(data=naturaldis2, ln_exports ~  deeppta + natdis_cont_deeppta | ID1 + ID2 + ID_dyad)
add_depth <-felm(data=naturaldis2, ln_exports ~  depth_index + natdis_additive | ID1 + ID2 + ID_dyad)

# Stargazer Table
stargazer(depth, enf, flex, contdepth, add_depth,
          header = FALSE,
          add.lines = list(c("Year Fe", "No", "No", "No", "No")), 
          omit.stat = c("adj.rsq", "f"), 
          model.numbers = FALSE,
          model.names = FALSE,
          dep.var.caption.include = FALSE,
          dep.var.labels.include = FALSE,
          column.labels = c("ln(Exports)", "ln(Exports)", "ln(Exports)", "ln(Exports)", "ln(Exports)"),
          title = "Natural disaster type and deep PTAs",
          notes = "",
          star.cutoffs = c(.05, .01, .001),
          covariate.labels = c("Deep PTA", 
                               "Natural disaster x Deep PTA", 
                               "Enforcement PTA", 
                               "Natural disaster x Enforcement",  
                               "Flexible PTA",
                               "Natural disaster X Flexible PTA",
                               "Cont natural disaster x Deep PTA",
                               "Additive Deep PTA",
                               "Natural disaster x Additive Deep PTA"
                               )
          )

```



### Table 3 ####
```{r xcaw}

# Creating interaction between deep pta and large disasters based on damage
naturaldis$natdislargedamage_deeppta <- NA
naturaldis$natdislargedamage_deeppta <- naturaldis$natdisdamage_large*naturaldis$deeppta
summary(naturaldis$natdislargedamage_deeppta)

# Interaction between PTA depth and large disasters based on damage
natlarge_reg <- felm(data=naturaldis, ln_exports ~  deeppta + natdislargedamage_deeppta | ID1 + ID2 + ID_dyad)



# Creating interaction between deep pta and medium disasters 
naturaldis$natdismeddamage_deeppta <- NA
naturaldis$natdismeddamage_deeppta <- naturaldis$natdisdamage_med*naturaldis$deeppta

# Interaction between PTA depth and medium disasters 
natmed_reg <- felm(data=naturaldis, ln_exports ~  deeppta + natdismeddamage_deeppta | ID1 + ID2 + ID_dyad)



# Creating interaction between deep pta and small disasters based on damage
naturaldis$natdissmalldamage_deeppta <- NA
naturaldis$natdissmalldamage_deeppta <- naturaldis$natdisdamage_small*naturaldis$deeppta

# Interaction between PTA depth and small disasters based on damage
natsmall_reg <- felm(data=naturaldis, ln_exports ~  deeppta + natdissmalldamage_deeppta | ID1 + ID2 + ID_dyad)



# Stargazer Table
stargazer(natsmall_reg, natmed_reg, natlarge_reg, 
          header = FALSE,
          add.lines = list(c("Fixed Effects", "Yes", "Yes", "Yes")),
          omit.stat = c("adj.rsq", "f"), 
          model.numbers = FALSE,
          model.names = FALSE,
          dep.var.caption.include = FALSE,
          dep.var.labels.include = FALSE,
          column.labels = c("ln(Exports)", "ln(Exports)", "ln(Exports)"),
          title = "Natural disaster size and PTAs",
          notes = "",
          star.cutoffs = c(.05, .01, .001),
          covariate.labels = c("Deep PTA", 
                               "Small natural disaster x Deep PTA", 
                               "Med natural disaster x Deep PTA", 
                               "Large natural disaster x Deep PTA")
          )


```


### Table 8 ###
```{r testing}

#Cont measure of natural disaster size (based on damage) interaction with depth 
lntotaldamage_deeppta <- naturaldis$lntotaldamage*naturaldis$deeppta

natcontinuous_reg <- felm(data=naturaldis, ln_exports ~  deeppta + lntotaldamage_deeppta | ID1 + ID2 + ID_dyad)



#Cont measure of natural disaster size (based on deaths) interaction with pdeth
lntotaldeath_deeppta <- naturaldis$lntotaldeaths*naturaldis$deeppta

natcontinuous_death <- felm(data=naturaldis, ln_exports ~  deeppta + lntotaldeath_deeppta | ID1 + ID2 + ID_dyad)



# Putting continuous measure of natural disaster size regressions together in stargazer for the Appendix
stargazer(natcontinuous_reg, natcontinuous_death, 
          header = FALSE,
          add.lines = list(c("Fixed Effects", "Yes", "Yes", "Yes")),
          omit.stat = c("adj.rsq", "f"), 
          model.numbers = FALSE,
          model.names = FALSE,
          dep.var.caption.include = FALSE,
          dep.var.labels.include = FALSE,
          column.labels = c("ln(Exports)", "ln(Exports)"),
          title = "Continuous natural disaster size and PTAs",
          notes = "",
          star.cutoffs = c(.05, .01, .001),
          covariate.labels = c("Deep PTA", 
                               "ln(Total Damage) x Deep PTA", 
                               "ln(Total Deaths) x Deep PTA")
          )


```


### Figure 1 ###
```{r setup, include=FALSE}

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

naturaldisfigure <- read_dta("naturaldisfigure.dta")


library(tidyr)
library(dplyr)

library(foreign)
library(haven)
library(ggplot2)
library(ggthemes) 


figure1 <- ggplot(naturaldisfigure, aes(x=year)) + 
  geom_line(aes(y=flood, color="Flood")) + 
  geom_line(aes(y=drought, color="Drought")) + 
  geom_line(aes(y=wildfire, color="Wildfire")) +
  geom_line(aes(y=earthquake, color="Earthquake")) + 
  geom_line(aes(y=extremetemperature, color="Extreme Temperature")) + 
  geom_line(aes(y=storm, color="Storm")) +
  geom_line(aes(y=epidemic, color="Epidemic")) +
  theme(panel.background = element_rect(fill = "white", colour = "white")) +
  theme_tufte() +
  scale_color_manual(values = c("black", "steelblue", "darkgreen", "blanchedalmond", "darkred", "navy", "darkgrey"), 
  name="Natural Disaster Type") +
  labs(x = "Year", y=  "Number of Natural Disasters") + 
  theme(axis.text=element_text(size=12),axis.title=element_text(size=12))

pdf("figure1.pdf", width=8, height=5.5) 
figure1
dev.off()

```



### Figure 4 ###
```{r}

nd_corr <- naturaldis %>% select(99:106) 
nd_corr_corr <- cor(nd_corr, use = "pairwise.complete.obs")
corrplot(nd_corr_corr, method = "color", type = "upper")

```


### Figure 5 ###
```{r}

# identifiying where the column is 
match("lntotaldamage",names(naturaldis))

damage_corr <- naturaldis %>% select(170:171) 
damage_corr_corr <- cor(damage_corr, use = "pairwise.complete.obs")

corrplot(damage_corr_corr, method = "color", type = "upper")

```


### Figure 6 ###
```{r dat, }

missdat <- read_dta("missingdata.dta")
cormat <- cor(missdat, use = "complete.obs")

cor_plot_miss <- ggcorrplot(cormat, hc.order = TRUE, 
                 type = "lower", outline.color = "white", 
                 lab = TRUE, lab_size = 2.5, tl.cex = 9,
                 ggtheme = ggplot2::theme_void)

ggsave("cor_plot_miss.png")

```





